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Abstract 

We check the universality properties of the two-dimensional Abelian sandpile model by 
computing some of its properties on the honeycomb lattice. Exact expressions for unit height 
correlation functions in presence of boundaries and for different boundary conditions are derived. 
Also, we study the statistics of the boundaries of avalanche waves by using the theory of SLE 
and suggest that these curves are conformally invariant and described by SLE 2 . 

PACS: 05.65+b, 89.75.Da 

1 Introduction 

Bak, Tang and Wiesenfeld introduced the theory of self-organized criticality as a general mechanism 
that can explain the behaviour of complex systems which naturally organize themselves into a 
critical state [1]. They denned the sandpile model as an example of slowly driven and dissipative 
complex system to explain the concept of self-organized criticality. Thanks to the Abelian property 
of the model [2], many statistical and dynamical results have been derived exactly. Among the 
main analytical results obtained for the isotropic two-dimensional model, one can mention 1-site 
probabilities of height variables [3], bulk correlations of height 1 variables and of some specific 
clusters known as weakly allowed clusters [4, 5, 6], the discussion of boundary conditions [7, 8] and 
the effect of boundaries on height probabilities [9, 10, 11, 12, 13], boundary correlations of height 
variables [14, 15, 16], bulk correlations of higher height variables [17] and avalanche and toppling 
wave distributions [18, 19, 20, 21]. The sandpile model has been investigated in the continuum limit 
and from a field theoretical perspective. It has been shown that the Abelian sandpile model can be 
described by a specific logarithmic conformal field theory, with central charge c = —2 [5, 22, 13, 17]. 

Most of the results, analytical or numerical, have been obtained when the sandpile model is 
formulated on the square lattice. In fact, the simple and symmetric structure of the square lattice 
make it easier to carry out some of the lattice calculations. However, other regular two-dimensional 
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lattices have been considered. Height probabilities and critical exponents of the sandpile model on 
the triangular and honeycomb lattices have been investigated using renormalization group trans- 
formations and numerical simulations [23, 24]. Moreover the critical exponents and the finite-size 
scaling functions for the avalanche wave distributions on the square, honeycomb, triangular, and 
random lattices have been evaluated from Monte Carlo simulations [25]. The results clearly suggest 
that the model on different lattices has the same set of critical exponents and scaling behaviour. 

In this paper, we study the Abelian sandpile model on the honeycomb lattice. The structure of 
the honeycomb lattice is more complicated than the square lattice, as there are two lattice points in 
each unit cell. However, the correspondence between recurrent configurations and spanning trees 
is maintained, and an exact expression for the Green function can be obtained, along with the 
exact asymptotic value. These are the two main ingredients to perform the lattice calculations and 
in turn, to check explicitly the universal behaviour of the model. We do that by looking at the 
universal terms of the unit height correlation functions with and without boundaries. 

The universality in the critical behaviour of sandpile model can be also checked from the point 
of view of geometrical features of the model. Indeed the dynamics of the model is such that each 
avalanche is formed on a compact domain with a boundary which converges to a fractal curve in the 
scaling limit [26]. For the sandpile model defined on the square lattice, it has been recently suggested 
that the boundary of these avalanche clusters belongs to the family of conformally invariant curves 
generated by the Schramm-Loewner evolution process SLE K , for a diffusitivity constant k = 2 [27]. 

Avalanches can be also decomposed into a sequence of simpler objects called toppling waves. 
While avalanches are believed to be described by a multifractal set of scaling exponents and have 
a complex scaling behaviour, waves show simple scaling properties and are more convenient for the 
analysis of avalanche statistics. Here, we also check the universality properties by studying the 
statistics of toppling wave boundaries for the model defined on the honeycomb lattice. 

This paper is organized as follows: in the next section we define the sandpile model on the 
honeycomb lattice and give the exact (and well-known) expression for the discrete Green function 
on the lattice. In the third section, we review the methods used for our lattice calculations and 
then give the details of the results in Section 4. We compare our results with the predictions of the 
c = — 2 conformal field theory in Section 5. Finally, we investigate in Section 6 the statistics of the 
toppling wave boundaries on the honeycomb lattice and verify its universality. We finish with some 
conclusions. The values of the lattice Green function for small distances are listed in an Appendix, 
which also contains the details for the calculations of its asymptotic behaviour for large distances. 

2 The sandpile model 

We start by defining the sandpile model on a two-dimensional honeycomb lattice of linear size L 
and with N sites. To each site of the lattice a random integer variable h G {1,2,3} is assigned 
which can be interpreted as the number of sand grains at that site. A configuration is characterized 
by the set of heights of all sites and is called stable if the height values are equal to 1, 2 or 3. 

The dynamics of the model is defined in two steps: 1) given a stable configuration, a grain of 
sand is added to a randomly chosen site, so that its height is increased by one, while the other sites 
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Figure 1: (a) A portion of the honeycomb lattice with the unit cells, marked by rectangles. Each 
unit cell contains two types of lattice sites, A and B. (b) The coordinate system used to label the 
cells is spanned by the two unit vectors a\ and 3,2- 

remain unchanged; 2) if the height of that site becomes greater than the critical height h c = 3, the 
site becomes unstable, topples and loses three grains of sand, each one of which drops on one of 
the nearest neighbours. The toppling rule can be written in the form hj — > hj — Ay for all sites j, 
where A« is called the toppling matrix and is equal to the discrete Laplacian, namely 



If, as a result of this toppling, some of the neighbours become unstable, the toppling process 
continues until all sites become stable and the avalanche ends. Thus in each time step, the dynamics 
takes the system from a stable configuration to another stable configuration. 

For a lattice with N sites, the total number of stable configurations is 3^, but not all of them 
occur in the steady state. The stable configurations are divided into two classes: recurrent and 
transient. After a long time, when the system enters the steady state, the transient configurations 
have zero probability of occurrence and all recurrent configurations occur with equal probability [2] . 
The burning algorithm allows to determine whether a given configuration is recurrent or not, and 
also establishes a one-to-one correspondence between recurrent states and spanning trees. Thus 
one can compute the total number of recurrent configurations by enumerating the spanning trees. 
From Kirchhoff 's theorem, the number of spanning trees is given by the determinant of the toppling 
matrix A, or the discrete Laplacian matrix. 

If we keep the form of the toppling matrix as above in (2.1), even for the boundary sites which 
have strictly less than three neighbours, then sand grains will leave the system each time a boundary 
site topples. Such sites are therefore dissipative and called open. If we set the diagonal element 
An equal to the number of neighbours of i, then i is conservative and called closed. Boundary sites 
can freely be chosen open or closed; the dynamics of the model, described above, is well-defined 
provided the system contains at least one open site. 




(2.1) 
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The honeycomb lattice can be divided into unit cells such that each cell contains two lattice 
points, one of type A (left site) and one of type B (right site), as shown in Fig. la. We choose the 
origin of the coordinate system at a site of type A. The position of a cell is specified by the vector 
na\ + ma2, where a\ and 0,2 are the unit vectors shown in Fig. lb and n,m are integer numbers in 
Z. Then, each site is characterized by the location of its cell, f= (n, m), and its location inside the 
cell, as a = 0, 1 for A-type or 5-type sites respectively. So the complete coordinates of a lattice 
site is a triplet (n,m;a). We will also use polar coordinates (n,m) <-> (r,<p) such that r = re tip , 
with r the Euclidean distance to the origin given by r 2 = n 2 + m 2 — nm and tp the angle measured 
counterclockwise from the h axis. Explicitly the change of coordinates reads n = r cos p + sin tp 
and m= ^= sin tp (so (r, ip) are the polar coordinates of the site (n, m; 0)). 

The diagonalization of A can be obtained by a method based on the decomposition of the 
lattice into unit cells [29]. One can write the connectivity between the vertices of the unit cells 
n = (ii,mi) and ti = (712,7712) in terms of 2 x 2 adjacency matrices a{f 1^2) given by 

,_ _ . f 1 if site «i in cell n and site 02 in cell r*2 are adjacent, .„ „. 

a aia2 (r 1 ;r2) = < " ' (2.2) 

I otherwise. 

As they only depend on the difference r = fz — r*i, we will write a(r) = o(ri;r2)- Explicitly the 
connectivity matrices of each unit cell with itself and its four neighboring cells read (rows and 
columns are labelled by a±, CK2 = 0, 1 in that order) 



«(0,0) = ^ Q J, a(l,0)=o(l,l)=a t (-l,0)=a t (-l,-l) = ^ 1 Q j, (2-3) 

while all the other matrices equal to 0. The Laplacian can then be written as 

A (fi;ai),(f 2 ;a 2 ) = [3 - a(0, 0)] <g> <S ft)?1 -0(1,0)®^^^ - o(-l, 0) ® 5 ft)?1 _ Sl 

- a(l, 1) (8) <V 2!f - 1+Sl+32 - a (-l, -1) (g) S^ ri _ Sl _g 2 . (2.4) 

Let us consider an 7V"i x N2 array of unit cells with periodic boundary conditions, namely the 
set of cells r = (n,m) with 1 < n < Ni, 1 < m < N2 and periodicity in both coordinates. In the 
decomposition (2.4), the matrices acting in the r space are made of cyclic matrices and are simul- 
taneously diagonalized by going to the basis of eigenfunctions Vfc l5 fc 2 (n,m) = e 2 ' mk ^ n / N ^ e 2nTk 2 m/N 2 _ 
In this basis, the Laplacian is block-diagonalized, A ~ (Be 1 ,e 2 M($i,02), with 

M(0i,0 2 ) = 31 -o(0,0) -a{l,0)e idl - o(-l, 0)e" iei - o(l, l)e i(9l+e2) - o(-l, -l) e -^+^) 

) , (2.5) 

where the angles are given by 9j = -j^-, with < k\ < N\ — 1 and < k2 < AT 2 — 1. 
From the previous results, we readily obtain 

det A = Yl det M (#i> ^2) = II [ 6 - 2 cos 6)1-2 cos 62 ~ 2 cos (01 + ^2)] , (2-6) 

01,02 01,02 
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a number clearly equal to zero since the eigenvector vo t o(n, m) = 1 is a zero mode. If however we 
leave out the zero eigenvalue, the Matrix- Tree theorem [28] states that the product over the non- 
zero eigenvalues, divided by the number of sites (= 2N1N2), yields the total number of spanning 
trees on the finite array of cells with exactly one site open. Since M(0, 0) has the eigenvalues and 
6, we obtain that this number equals 

— — — 6-2cos — 2 cos— 2 cos — h . (2.7) 

(fci,fe 2 )^(0,0) 

From this, one obtains the entropy per site for the sandpile model on the honeycomb lattice in the 
thermodynamic limit [29], 

7T J/3 J/1 

1 9 2 log 6 - 2 cos 6>i - 2 cos 6> 2 - 2 cos (6>i + 6> 2 ) - 0.807665. (2.8) 

87T Z L J 

The effective number of degrees of freedom per site in the set of recurrent configurations is thus 
equal to e z ~ 2.243. 

On the infinite lattice, the Laplacian (2.4) depends on r = r 2 —r\ only, A Q1)Q , 2 (r) = A^.q,^^.^), 
and reads in Fourier space, 

A ai , a2 (6) = £ A Q1 , Q2 (0 e^= 31 - 5> Ql , aa (r) = M^^, 2 ), (2.9) 
f r 

where now = (9\,6 2 ) belongs to [0, 2ir] 2 . The Green function is thus the inverse Fourier transform 
of the inverse of M and also depends on r only, 

d6i d9 2 1 



f7r cl^i d^2 e —i(n2—ni)0i e —i(m2—mi)d2 



n 4vr 2 6 -2 cos X - 2 cos 2 - 2 cos (6>i + 2 ) 

/ 3 i + e -*i +e -i(*i+fc)\ 

X ^i + e^+^+fe) 3 J' ( 2 - 10 ) 

From this it immediately follows that Gaa(0 = GBB{r) and G A b(t) = Gea^—^), as well as 

o /•/•7T e in9i+im02 

G AA (n,m) = -^ d9 1 d6 2 - (2.11) 

87H jj— 71- 3 — cos pi — cos 6> 2 — cos (tq + 6> 2 ) 

G AB (n,m) = -[G AA (n,m) + G AA (n+l,m) + G AA (n + l,m+l)}, (2.12) 

G B A(n,m) = -[G AA (n,m) + G AA (n-l,m) + G AA (n-l,m-l)], (2.13) 

in agreement with Poisson's equation. 

For the lattice calculations developed in the next sections, we need to know the Green function 
for small distances, and its asymptotic behaviour for large distances. These are discussed and 
collected in the Appendix. 
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3 Lattice calculations : general formalism 



Our main purpose in this paper is to check the universality properties of the sandpile model by 
computing some of its properties on the honeycomb lattice, for instance joint probabilities of local 
height variables. Among these the simplest are related to so-called weakly allowed clusters (WACs). 
These are finite clusters of sites with specific height values such that decreasing the height of any 
of their sites by one turns them into forbidden subconfigurations [2, 4]. 

The method used to compute the joint probability of such clusters, due to Majumdar and Dhar, 
is based on a modification of the toppling matrix [4] . They showed that the number of recurrent 
configurations which contain a given WAC is equal to the total number of recurrent configurations 
of a new sandpile model, defined in terms of a modified toppling matrix. The modification is 
usually done by removing some connections from the cluster to the rest of the lattice, and at the 
same time, by adjusting the diagonal elements of the toppling matrix in order to prevent the sites 
from being dissipative. Then the new toppling matrix can be written in terms of the original one 
as A new = A + B, where the defect matrix B encodes the modification: Bij = Bji = 1 if the 
symmetric bond between sites i and j is removed, and Ba = — n, if n bonds have been cut from the 
site i. Since the modifications concern a finite collection of sites, the matrix B has finite rank. The 
probability of occurrence of a weakly allowed cluster S is obtained by computing the determinant 
[4] 

Hpt Anew 

P(5) = ^eI^ = det(l + GjB) ' (3 - 1} 

where G = A -1 is the lattice Green function. Because B has finite rank, the matrix indices in the 
determinant may be reduced to the sites affected by the modification, so that the determinant is 
actually finite. 

The same method can be used to calculate the probability of occurrence of several WACs, and 
therefore their correlation functions. Each cluster Si comes with its own defect matrix Bi, and 
the full defect matrix is simply the direct sum of all B^s. The matrices B and G acquire a block 
structure, where each block refers to the sites involved in each of the clusters. For example, a 
2-cluster correlation function can be found by computing 

g( B ; D). ,3.2, 

If Si and S2 are respectively located around r\ and r*2, this probability will depend on Green 
function entries for small distances within S\ and Si (those entries in G\\ and G22), and on 
entries labelled by pairs of sites, one being close to n, the other being close to r*2. As one is mainly 
interested in correlations of clusters separated by large distances, the calculation of the determinant 
requires to know the Green function for both small and large distances. Calculations of multicluster 
probabilities, for about a dozen of different WACs and for up to four clusters, have been carried 
out on the square lattice [5]. 

Note that the above formalism remains valid for a finite or infinite grid, with or without bound- 
aries. The form of B will in general depend on whether some of the clusters touch a boundary; in 
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addition the Green function used in the calculations must be appropriate to the geometry and the 
boundary conditions used. The thermodynamic limit can be conveniently evaluated by using the 
Green function on the infinite lattice directly. 

The simplest of all WACs consists of a single site with height equal to 1. In this case, the 
modification consists in leaving only one bond between the height 1 and the rest of the lattice, and 
correspondingly by decreasing the diagonal entries of A as explained above. On the honeycomb 
lattice, it means that B is identically zero except for a 3-by-3 block 



labelled by the site where the height is fixed to 1, and any two of its nearest neighbours. If the site 
where the height is 1 is on a boundary, the corresponding non-zero block is smaller. 

^From correlations of the above kind, one may infer the scaling behaviour of lattice observables 
like the height 1 variable in the bulk or on a boundary with a given boundary condition. In the 
scaling limit, the correlations in the bulk, on a boundary, or at a finite but large distance from 
a boundary, are all universal quantities, which the underlying conformal field theory is able to 
describe. Being universal, they should not depend on the type of lattice on which the model is 
defined. This is what we want to check in the following. 

On the square lattice, it has been shown [4] that the height 1 variable in the bulk scales like a 
dimension 2 field, implying in particular that the 2-point correlation decays like 1/r 4 . The same 
is true for the height 1 variable on an open or a closed boundary [14, 15]. Once the scaling limit 
of the height 1 variable has been properly identified with a specific field of the conformal theory, 
higher correlation functions are fixed. On the square lattice, the lattice results confirm well the 
predictions of logarithmic conformal field theory [12]. 

In contrast, higher height variables in the bulk do not scale the same way since their correlations 
involve logarithmic functions of the distances [12, 13, 17]. For instance the 2-point correlation of a 
height 1 variable and a height variable strictly bigger than 1 has been shown to decay like log r/r 4 
[17], whereas that of height variables greater or equal to 2 is conjectured to decay like log 2 r/r 4 [13]. 
On an open boundary, all four height variables scale the same way, while on a closed boundary, 
the height 2 and 3 variables have a slightly different but still non-logarithmic scaling compared to 
the height 1. All known results on the square lattice are consistent with the identification of the 
height 2, 3 or 4 variables in the bulk as the logarithmic partner of the height 1. 

Likewise on the honeycomb lattice, height 2 and 3 variables are expected to have a logarithmic 
scaling similar to the heights 2, 3 and 4 on the square lattice. Here we will restrict ourselves to the 
height 1, and check that it can be identified with the same conformal fields as on the square lattice. 

4 Lattice calculations : results 

The calculation of height 1 probability is the simplest case. We give in this Section the results 
we have obtained for the various correlations, in the bulk with and without boundaries, and along 
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boundaries. 
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4.1 In the bulk 



As recalled above, if we want to have a height 1 at a site i, we should remove the bonds between i 
and two of its neighbours, and change the diagonal entries of the toppling matrix. This corresponds 
to take the defect matrix B as above in (3.3). This defect matrix is non-zero on three sites only, 
namely i and the two chosen neighbours, so that the general formula requires to know the Green 
function on the same three sites. From the results of the Appendix and obvious symmetries, it 
reads, for the infinite lattice (we assume i is of type A), 

fG AA (0,0) G AB (0,0) G AB (0,0)\ ( ~\ ~\\ 

G = 



3 3 

1 o — - 

3 u 2 



G AB (0,0) G aa (0,0) GaaO-,0) I = G AA (0,0) + | -i 

3 ~2 



(4.1) 



\G AB (0,0) G AA (l,0) G AA (0,0)/ \-\ ~h i 



The constant piece G AA (0, 0) is divergent on the infinite lattice, but drops out in the product GB 
since B has column (and row) sums equal to 0. We thus find that the probability that a given site 
has a height equal to 1 is given, in the infinite volume limit, by 

Pi = det(I + GB) = j2~ 0.0833. (4.2) 

The other two single height probabilities are known from numerical simulations [25] , and found to 
be P 2 ~ 0.2937 and P 3 ~ 0.623. 

We can similarly compute multi-site height 1 probabilities. We first consider the two-site prob- 
ability for two A- type sites, one at the origin and the other at (n, m;0). The two-site height 1 
probability is obtained according to the formula in Eq. (3.2), for which the Green function for sites 
close to the two reference points is required. This is computed for large distances and arbitrary 
positions in the Appendix, and we obtain 1 

„ / ^ „9 /- 3 1 1 2 + 5cos6<p \ . , . 

Pn(n,m; 0) = P? (l - - — 2 - & +...). (4.3) 

The two-site probability when the distant site is of type B can be computed in a similar way, 
and reads 

^ , \ „9 / 31 12 — 5 cos 6(fi \ , , , . 

Pn(n, m; 1) = P^l - - _ —JL + . . . j , (4.4) 

where r and ip are the polar coordinates of the site (n, m; 1). We see that the subdominant term 
~ r~ 6 not only depends on the angular position of the distant site but also on its type, A or B. 
This is a first sign that only the dominant term in r~ 4 is universal, and rotationally invariant, as 
expected for a scalar observable. 

In both cases, the 2-point correlation of two heights 1 separated by a large distance r behaves 

like 

P llW -Pf = + (4.5) 



1 The analogous calculations on the square lattice [5] have been carried for specific spatial configurations of the 
heights f , mostly when they are aligned on a principal or diagonal axis. The use of the Green function for arbi- 
trary positions from the Appendix enables us to read off more clearly the universal terms and allows a more direct 
comparison with conformal field formulas. 
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Up to the numerical coefficient, this is the same result as on the square lattice: interpreted, in the 
scaling limit, as the 2-point correlation function of the field <f> associated to the presence of a height 
1, it implies that is a scalar field with scaling dimension 2, and fixes its normalization. 

The general three-site probability for three heights 1 can also be computed. Here we restrict 
ourselves to ^4-sites only, in cells located at r\ = 0, r*2 and The probability depends on the 
vectors fij = r. \ — fj which we write in polar coordinates as f*y = r^e 1 ^ , so that ipij = it + tpji. 
We find the following result for the connected 2 three-site probability, 

„ ,-. -. 1 1 ( sin 2(9912 - v?i 3 ) cos 3<£>23 

All(ri,r2,r3) CO nn =~ 



5767T 3 (ri 2 ri 3 r2 3 ) 2 I r 23 

+ sin 2(1^12 - y? 23 ) cos3y9i 3 + sin 2((p 13 - (^3) cos 3y 12 1 ^ g . 

ri3 r i2 j ... . 

The main observation is the absence of a dimension 6 term, which, in the scaling limit, should 
correspond to the 3-point function of the field 0. Indeed one sees that on the lattice, the dominant 
term of the three-site correlation has scale dimension 7 (and even 8 for specific spatial configura- 
tions) . The same observation was made on the square lattice [5] , and leads to the expectation that 
the 3-point function of (f> vanishes identically. 

Finally, in the same notations, we have computed the general 4-site probability, again for A-sites 
only. The connected part reads 

D / -» -» -» \ 1 /\/3-Pl\ 4 [COS 2(9313 - (f 14 - ip 2 3 + V?24) 

-pLlll(ri,r2,r3,r 4 )conn = -t( <^ -, 7n 

4 V 7T / I (ri 3 ri4r2 3 r 2 4) 

COS 2(93i2 - 9314 ~ V?23 + ^34) | COS 2(y3i 2 - V313 - y3 2 4 + V?34) \ | ^ ^ 

(ri2n4r23^34) 2 (n2n3^24r34) 2 / 

As on the square lattice, the dominant term has the expected scale dimension 8, and should 
correspond to the 4-point correlator of (f>. 

4.2 On upper-half planes 

In addition to probabilities on the infinite lattice, height 1 probabilities on semi-infinite lattices can 
also be examined. We will consider two upper-half planes, one bordered by a tilted boundary, the 
other by a horizontal boundary, as shown in Fig. 2 and Fig. 3. The tilted boundary is parallel to the 
h axis, and, for this reason, will be called principal; it is made of all A-sites on the line m = 1; each 
boundary site has two neighbours which lie in the interior of the upper-half lattice. The horizontal 
boundary contains the sites in all the cells on the line n = 2m — 2; each boundary site has again two 
neighbours, one of which is itself a boundary site. On each type of boundary, the uniform open or 
closed condition corresponds to set An = 3 or Ajj = 2 respectively for all boundary sites. Although 
the technical details differ for the two kinds of boundaries, the results should not, and the joint 
probabilities should only depend on the distances separating the heights 1 and the boundary. 



2 The connected n-site probability is obtained by subtracting from the full n-site probability products of lower 
order probabilities. 



9 



(a) (b) 
Figure 2: Principal boundary with closed (a) and open (b) boundary condition. 



The calculation of height 1 probabilities on a upper-half plane follows the same principles as on 
the infinite lattice. In particular the same defect matrix B can be used if the height 1 is not at a 
boundary site. The only difference is that we should use the discrete Green function adapted to 
the boundary condition we choose. Open and closed Green functions are usually obtained using 
the image method, and, except in one case, the same method may be used here too. Let us first 
consider the principal boundary, in Fig. 2. 

For the closed boundary condition, each site of the half-plane above or on the boundary is 
mirrored through the reflection line shown as the dotted red line. An A-site has a mirror image 
which is a B-site and vice-versa. The coordinates of the images are indicated in Fig. 2a. The closed 
Green function for the principal boundary then reads 

^(ni,mi;a)(ii2,m2;0) — ,m\ ;a)(n 2 ,m 2 ;0) + G{n\ ,m\ ;a)(n2—ra.2, 1—7712 ;1) ' (4-8) 

^(ni,mi;a)(n2,m 2 ;l) — ^ {n\,mx-,<x){n%,m%V) ^(ni,mi;a)(n2— ma+1,1— W2;0) ■ (4-9) 

For the open boundary condition, the situation is slightly more complicated. The reflection line, 
shown in Fig. 2b as the dotted red line, is such that a B-site above the boundary reflects itself to a 
B-site below the boundary, however the mirror image of an A-site does not belong to the lattice. 
Instead the strict mirror of an A-site is the center of an hexagon lying below the boundary. We 
then define the three B-sites on this hexagon as the three mirror images of the A-site above the 
boundary, each one being weighted by a factor 1/3. The open Green function is then related to 
the Green function on the plane by the following expressions, 

1 



LT (ni,mi;a)(n 2 ,m 2 ;0) ^(ni,mi;a)(n 2 ,m 2 ;0) ^ 



^(ni,mi;a)(ii2-m2,-m2;l) 
^(ni,mi;tt)(«2- m 2 — 1,— m 2 ;l) ^'(ni,mi;o)(n 2 — m 2 ,— m 2 +l;l) > (^-10) 
^(ni,mi;a)(n 2 ,m 2 ;l) = ^(ni ,mi ;o)(n 2 ,m 2 ;l) — ^(ni ,mi;»)(n 2 -m 2 ,-m 2 ;l) • (4-11) 
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Figure 3: Horizontal boundary with the open boundary condition. 



These formulas and the asymptotic behaviour of the bulk Green function enable us to obtain 
the asymptotic behaviour of G op and G cl , and in turn, the scaling form of the height 1 probabilities. 
The height 1 probability at an arbitrary point (n, m; a) does not depend on n, by translational 
invariance along the h axis. When m S> 1 and for a = 0, we obtain the following 1-site probabilities, 

P°P(m) -P 1 = +—^\ + —^=- \ + ... = A + • • • . ( 412 ) 

18\/3vrm 3 vr 4r 2 

r i . . 11 11 \/3Pi 1 . . 

lf{m)-P 1 = = 5 = _ + ... = _* L +... 4.13 

where r = \/3(m — l)/2 is the Euclidean distance from the site (n, m;0) to the boundary. The 
distinctive change of sign between the two types of boundary conditions was also found on the 
square lattice [9, 5]. 

Let us now look at the second, horizontal boundary, shown in Fig. 3. We try to adapt the image 
method to this boundary in order to compute the discrete Green function for open and closed 
boundary. For the open condition, the image points are found by a reflection with respect to the 
horizontal line n = 2m — 1, pictured as the dotted red line in Fig. 3. Such a reflection preserves the 
type, A or B, of sites, so that we obtain 

^(ni,mi;a)(n 2 ,m 2 ;/3) = ^(ni,mi;a)(n 2 ,m 2 ;/3) — G ! (ni,mi;a)(n2,ri2-m2+l;/3) ■ (4-14) 

In contrast, for the closed boundary condition, we have not been able to define the appropriate 
reflection map, with one mirror image or several mirror images like above, so as to make the image 
method work 3 . As a consequence, we could compute the 1-site height 1 probability for the open 



incidentally, the diagonal boundary on the square lattice, defined by the line x = y in Z 2 , may be considered in 
addition to the more usual boundaries parallel to the principal axes. In this case however, the image method works 
for both the open and closed boundary conditions. To our knowledge, no calculation has been carried out with this 
type of boundary. 
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condition only. The probability for finding a height 1 at a site [2y — 2, y + m — 1; a), located at a 
distance r = m — 1 from the horizontal boundary, does not depend on y. So we may put y = 1 and 
consider the site (0,m;a), in which case the probability reads 

pr p (m) - p l = +—^\ + ... = A + • • • ( 4 - 15 ) 

with the same dominant term as for the other boundary, as expected. 
4.3 On boundaries 

Finally the height 1 correlation functions along a boundary can be computed for the types of 
boundaries considered above. In each case, a boundary site has only two neighbours so that the 
defect matrix required to force a height 1 is two-by-two and depends on the boundary condition, 
open or closed. Explicitly they read 

On the principal boundary, it is not difficult, with the appropriate Green functions given in the 
previous subsection, to compute the one-site probabilities for a height 1, on an open or a closed 
boundary, 

F? = H + -±- - -I ~ 0.129, Pf = ^ - \ ~ 0.218. (4.17) 
36 V3vr 7r 2 7r 3 

We note that the height of a closed site can only take the values 1 and 2, so that the previous result 

implies P 2 cl = § - & ~ 0.782. 

For the joint probabilities of having two, three or four heights 1 separated by distances along 

the principal boundary, we obtain, for the open boundary condition, 



Pi? 



Pill 
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+ . . . (4.20) 
and similar expressions for the closed boundary condition, 
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= _I(^fi_ + ... (4.21) 
4 V 2vr / rf 2 V J 

l/\/3\3 1 

= _( ) |_ _ _ _ (4.22) 
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On the horizontal boundary, and for the open boundary condition only (as explained above, we 
cannot handle the closed condition) , we obtain similar results for the single site height 1 probability, 

oy o q 

and for the first correlation functions, 



= -j(- 2 -f) \ 1 , 2 +--- (4-26) 

conn 4 \7T Z 97T / \r12r13r23) 

poP _ i/i An 1 1 1 1 1 ] 

1111 conn 8V7T 2 9fT > \ {r^Tl^iTu) 2 (ri2^14^23^34) 2 (nsrU^S^) 2 J 

+ . . . (4.27) 



5 Conformal Field Theory 

If the above results for heights 1 on the honeycomb lattice are compared with those obtained on 
the square lattice [5, 6, 15, 16], it is immediately clear that they coincide, up to normalizations, and 
thereby confirm the universality of the field assignements. So we restrict here to a brief reminder 
of the main features of the conformal field interpretation of these lattice results, and take the 
opportunity to collect the various formulas. 

On the square lattice, it has been shown that, in the scaling limit, i.e. in the large distance 
limit, the joint height 1 probabilities on the lattice are exactly reproduced by conformal correlators 
of primary fields. If the height 1 variable lives in the bulk, the corresponding primary field is 
a non-chiral field with conformal weights (1,1), whereas it is a chiral, boundary primary field of 
weight 2 in the case the height 1 lies on a boundary, closed or open. It turns out that all these 
correlators can be understood, and computed, by writing the primary fields in terms of a pair of 
symplectic free fermions 6,9. 

The theory of symplectic free fermions, with central charge c = —2, is the logarithmic conformal 
field theory which is, by far, the best understood, see for instance [30], and the more recent work 
[31] as well as the references therein. We only need here the most basic features of it. 

The symplectic fermions are anticommuting, space-time scalar fields with propagators given by 

9(z,z)d (w,w) = 6(z,z)6 (w,w) = 0, (5.1) 
9(z,z)9(w,w) = — \og\z — w\, (5.2) 

from which all higher correlators may be obtained from Wick's theorem. Since the propagator 
(5.2) is a sum of a chiral term and a antichiral term, the fermions satisfy dd6 = dd6 = in 
all correlators. A Langrangean realization is provided by the action S ~ J dOdO, which has the 
previous two conditions as equations of motion. 

The fermion fields themselves are not primary, but their first derivatives are primary. In par- 
ticular (j)(z,z) = dd(66) is a primary field with conformal weights (1,1). It is not difficult to 
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compute its 2-, 3- and 4-point functions. They are given explicitly by the following expressions 
where zij = Zj — Zj , 



(mm) = (5-3) 



i 

2M 

(0(1)0(2)0(3)) = 0, (5.4) 
(0(1)0(2)0(3)0(4)) = + ^ + ^ 

If 1 | 1 | 1 \ 

8 I (Zl2^34^13^24) 2 (^13^24^14^23) 2 (^14^23^12^34) 2 ' 

The 3-correlator vanishes identically because the various Wick contractions necessarily involve the 
contraction of 89 with 89, or 86 with 89. In the 4-point correlator, the first three terms are products 
of 2-point functions and are not part of the connected correlator. 

The correlation functions of on the upper-half plane can be similarly computed by using the 
appropriate Green function, namely 

9(z, z)6(w, w) = — log \z — w\ ± log \z — w\, (5-6) 

with the + sign for the closed boundary, and the — sign for the open boundary. It yields in 
particular the 1-point function of on the upper-half plane, 

1 

"4(Im2:) s 



A simple comparison with the results obtained in the previous section shows that the leading 
terms of the connected joint probabilities are exactly reproduced by the above correlators provided 
the subtracted height 1 variable 5 (hi — 1) — P\ converges, in the scaling limit, to acf)(z, z) for some 
normalization a. The results in the bulk imply that for the honeycomb lattice, a = ± ■ The 
results on the upper-half planes then fix the sign, 

q , c . _ _ vmi _ i m 



vr 4^- 



7T 

By comparison, the results for the square lattice imply a sq = — P{* q = — 2 n a [5]. The way this 
specific conformal field emerges in the scaling limit has been demonstrated in [22]. Moreover, from 
the conformal field theory point of view, the open and closed boundary conditions have been shown 
[7] to be related to each other by the insertion of a chiral primary field of conformal weight —1/8, 
and leads to the change of sign in the 1-point function of 4> on the upper-half plane. 

For the purpose of describing the boundary height 1, we need the chiral version of the previous 
fields. So one also considers chiral symplectic free fermions with contractions 

9(z)9(w) = 9(z)9(w) = 0, (5.9) 

9(z)9(w) = ~\og(z-w). (5.10) 
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The chiral version of <f> that we will use, namely (f> c = d9d9, is not a primary field since it is 
proportional to the stress-energy tensor of the Lagrangean realization, T(z) = 2d8d6. The first 
correlators of <p c with itself read 

<&(1)&(2)> = - jV' ( 5 - n ) 

42 12 

{M1)M 2 )M3))= -i J -l— ? , (5.12) 
(^(1)^(2)^(3)^(4)) . + + 

~a{? _ ^2+7 " ^2+7 " ^1- (5 ' 13) 

8 1(2)12^13^24^34) {Z12ZUZ23Z34) (^13^14^23^24) J 

The boundary 2-, 3- and 4-correlators computed in Section 4.3 have exactly these functional 
forms, and show that the boundary height 1 variable, subtracted with the appropriate value of P±, 
converges to a c (/) c . The normalization depends on the type of boundary, principal or horizontal, 
and on the boundary condition. One finds 



a h.c.,princ,op = L a h.c.,prmc,cl = _ 

2y/3ir vt 2 2tt 



^i.c.pvijicop _ jj; ^ ^h.c.jprinc^l 

2^/3^ 7 

h.c.,horiz,op = VO 

vr 2 9vr v ; 

The two normalization factors for the open condition are positive, whereas the normalization for 
the closed condition is negative. 

On the square lattice, the boundary height 1 variable was also seen to converge to 4> c with a 
normalization, on a boundary parallel to a principal axis, given by [15, 16] 

«— l-iS*^ <**--§(!-§)■ 



Again the normalization is positive for the open, and negative for the closed boundary condition. 

It should be emphasized that, whereas the scaling limit of the height 1 variables, in the bulk 
and on open/closed boundaries, can be described by conformal fields which are themselves related 
in a simple way to symplectic free fermions, it is not so for all observables of the sandpile model. 

On the square lattice, it has been shown that boundary higher height variables scale to conformal 
fields which have simple expressions in terms of symplectic fermions. On an open boundary, the 
heights 2, 3 and 4 scale to the same field (f) c as the height 1, while the heights 2 and 3 on a closed 
boundary have a slightly different scaling, since they converge to a combination of <p c and 9dd9. On 
the honeycomb lattice, only the field 4> c is expected (a closed site has two neighbours, and therefore 
its height takes only two values). 

In contrast, the higher height variables in the bulk are all described, up to normalization, by 
a single scaling field ip, which turns out to be a logarithmic partner of the field 4> describing the 
height 1 in the bulk. However the reducible but indecomposable representation they generate does 
not belong to the theory of symplectic fermions 4 . Whether this representation has a Lagrangean 



4 The recent article [32] is a general study, in a much broader context, of classes of (chiral) representations such as 
the representation generated by the pair 4>,tp, which appears in their Example 7. 
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realization is an open and important problem. The same distinction between the height 1 and the 
higher heights in the bulk is expected on the honeycomb lattice, or indeed on any regular lattice. 



6 Boundary of Wave Clusters and Conformal Invariance 

The scaling behaviour of the two-dimensional critical lattice models can be reflected in the statistics 
of non-crossing random curves which form the boundaries of clusters on the lattice. In the 1920's, 
Loewner studied simple curves growing from the origin into the upper-half plane EI [33]. Loewner's 
idea was to describe the evolution of these curves in terms of the evolution of the analytic function 
g t , which maps conformally the region outside of the curve into H. He showed that this function 
satisfies the following differential equation 

for a real continuous function £t, related to the image of the tip of the curve under gt. Conversely, 
a continuous real function £ t implicitly defines a curve growing in H. 

Much more recently, Schramm followed the idea that a measure on the continuous driving 
functions ^ would induce a measure on the set of growing curves in H, and showed that the latter 
measure is conformally invariant if and only if the former measure is the Wiener measure for the 
standard one-dimensional Brownian motion Bt [34]. This subsequently led to a completely new 
perspective on random curves arising in conformally invariant critical systems, see [35] for a review. 
In this context, setting = \JHBt for different parameter k corresponds to different universality 
classes of critical behaviour. 

Avalanche boundaries in Abelian sandpile model are random dynamical curves whose statistics 
can be studied using the theory of SLE [27]. It has been suggested, on the basis of numerical 
simulations on the square lattice, that the boundaries of avalanche clusters are conformally invariant 
with the same properties as loop erased random walk model, with diffusivity constant k = 2. 

Since an avalanche has a complicated structure, understanding its dynamics can be simplified 
by decomposing the avalanche into a sequence of more elementary objects called toppling waves 
[18]. The waves are constructed as follows. If, as a result of the addition of a grain to a site i, that 
site i becomes unstable, it topples, as do the sites which become unstable as a consequence of the 
first toppling at i. The first wave is the collection of all sites which have toppled given that the 
initial site is not allowed to topple more than once. One can show that the sites in the first wave all 
topple exactly once. After the first wave is completed, the initial site, if still unstable, is allowed to 
topple a second time, and doing so, triggers the second wave of topplings. The process continues, 
with a third wave, fourth wave and so on, until the initial site i is stable and the avalanche stops. 
The important property of waves is that they are individually compact (no hole), and the sites in 
each wave topple exactly once. 

To check the universality of the model, we consider an ensemble of wave boundaries, on the 
honeycomb lattice, and repeat the analysis carried out in [27] for the square lattice. At first, we 
calculate the fractal dimension df for the wave boundaries, determined by the scaling relation 
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Figure 4: Main frame: Log-log plot of the length of wave boundaries I versus the radius of gyration 
R, simulated on the honeycomb lattice with the linear size of 2048. Inset: Log-log plot of the 
average area of wave clusters A versus the length I. 



I ~ R d f between the perimeter I of the curve and the radius of gyration R. The result for waves is 
df = 1.25 ± 0.01, see Fig. 4. The inset of Fig. 4 shows the scaling of the mean area of the wave 
clusters with their perimeter length as A ~ l 2 / d f , which is consistent with the one discussed in [36]. 
Furthermore, from the relation df = 1 + § for the fractal dimension of SLE curves [37], this fractal 
dimension is consistent with the value k = 2, obtained for the boundary of avalanche clusters on 
the square lattice [27]. The central charge associated with k = 2 is c = (3h,-8)(6-K) _ _2 jggj_ 

One of the questions about SLE curves that has a neat answer is the following: for a curve 
connecting two points on the boundary of a domain D, what is the probability that the curve 
passes to the left of a given point interior to the domain ? It is usual to take the domain D to be 
the upper-half plane and the boundary points to be the origin and the point at infinity. In this 
case, an interior point of the domain is represented in polar coordinates as z = Re 1 ^ . By scale 
invariance, the above probability depends only on 4> S [0, tt] and is given by [39] 

pm = \ + 7^%y ^ Q; % - ^ ( 6 - 2 ) 

where F12 is the hypergeometric function. For k = 2, this reduces to 

20 - sin 2(f) 

P 2 (0) = 1 ^ • (6.3) 

In order to check Eq.(6.3) for wave curves (loop curves), at first step we should convert these 
curves to curves which connect the origin to the infinity (chordal SLE). To this aim, we cross any 
given loop by an arbitrary straight line as real line at two points xo = and and consider only 
a segment of curve which is above the real line. Then by the conformal map (p(z) = *°°* z , curves 
in the upper half plane are transformed to a set of curves connecting the origin to infinity. 
The computed probabilities for points at distances R = 0.1,0.2,0.4 and 0.5 is consistent with 
Eq.(6.2), with k = 2.1 ± 0.1 (see Fig. 5 (a)). 

A more careful test which shows the correspondence with SLE, is to extract the Loewner driving 
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Figure 5: (a) The probability that a wave boundary passes to the left of a point with polar 
coordinates (R,<p), for R = 0.1,0.2,0.4 and 0.5. The solid line shows the function P2(<ft) obtained 
from SLE for k = 2. (b) Statistics of the driving function for the wave boundaries in the 
sandpile model on the honeycomb lattice. Main frame: the linear behaviour of (£(i) 2 ) with the 
slope k = 2.0 ± 0.2. Upper-right inset: the diffusion coefficient is k = 2.0 ± 0.2. Lower-left 
inset: the probability distribution function of the noise different colors correspond to 

t = 0.02,0.04,0.06,0.08. 



function There is an algorithm for chordal SLE curves, based on the approximation that driving 
function is a piecewise constant function [40]. As we mentioned at previous case, with conformal 
map ip{z) = z°°l z i the- segment of loop curves at upper half plane are converted to chordal curves. 
Then, each curve is parameterized by the dimensionless parameter t (that it is not time). In this 
case the Lowener equation is as dg t /dt = 2/{ip' (g t )[ip(g t ) — £ t ]}, which gt{z) maps the half-plane 
minus the curve up to t into the whole of upper half-plane. For a constant £, the equation can be 
solved for g t as: 



n f ^_„ Wooixoc - z) + yggz - i])' 2 + 4t(xoo - z) 2 x (xqq - r?) 2 

Lr t,£\ z ) — x oo — 7T-, : 1 a i \o 777 no / no 

^Sol^oo - Z) + V x 5o{ Z ~ V) + 4t(Xoo - Z) 2 X (Xoo - T)) 1 

Where, r\ = y? -1 ^). 

According to algorithm, the interval [0, T] is divided to smaller intervals [t n ,t n +i) with to = and 
^7V+i = T, such that the is approximated by the constant ^ n = £(t n ) in each interval. In this 
case, the function of gt is expressed as composition of G tN ~t N _ 1 £ N _ 1 o ■ ■ ■ oG tl £ . The action of each 
of Gt^(z) is such that when they apply on the points of the curve, remove the first point from the 
sequence of points. 

Now, we follow this algorithm for extracting of the driving function. At First, we take the 
points of our curves on upper half plane approximated with {zq = 0, z±, . . . , Zn = a?oo}- Then, 
using the parameters rjo = V 9_1 ('?o) = [ReziXoo — {Rez\) 2 — (Imzi) 2 ]/(x 00 — Rez\) and t\ = 
(ImziYx^/ '{4[(Rezi — x^) 2 + (imzi) 2 ] 2 },the map G tl ^ applied to the points resulting in a new 
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sequence, by one element shorter: z' k = £ (zfc+i), with k = 1, . . . ,N. The operation is iterated 
on the new subsequence of points until one obtains the full set of and for each curve. The 
result of this procedure is an ensemble of £(t) that its statistics as shown in Fig. 5 (b) , converges 
to a Gaussian process with variance (£ 2 (i)) = nt and k = 2.0 ± 0.2. This result, together with the 
other evidences, certify that the wave boundary curves of sandpile model are conformally invariant 
and described by the SLE2. 

This result seems to be reasonable from the correspondence with the spanning trees. In fact, one 
can construct exactly a two-component tree on the lattice, representing a wave [18]. The bound- 
ary of a wave as the dual of the spanning tree is expected to be SLE2. The statistics of the wave 
boundaries with diffusion coefficient k = 2 confirm the relation of the sandpile model with a c = — 2 
conformal field theory. 



7 Conclusions 

In this paper, we have investigated some of the properties of the Abelian sandpile model on the 
honeycomb lattice. The scaling behaviours of the height correlation functions in the bulk, in the 
presence of boundaries, and on boundaries, are in the agreement with those obtained on the square 
lattice, and correctly predicted by a c = — 2 conformal field theory. 

We have also checked the universality properties of the model from the point of view of its 
geometrical features, namely the statistics of the boundaries of the toppling waves. We found 
numerically that the boundaries of wave clusters are conformally invariant, and well described by 
the SLE process with diffusivity k = 2. 



A Appendix 

We collect in this Appendix some of the values of the lattice Green function on the honeycomb 
lattice, for small distances, and also give its asymptotic behaviour for large distances. 

In the coordinate system used in Section 2, the Green function on the honeycomb lattice for a 
pair of points of the A type and separated by the vector (n, m), is given by 

G AA (n,m) = — // dflidto ^ 7. 2 7a T~2j r (A.l) 

8ir z JJ-TT 3 — cos 0i - cos 6 2 — cos (0i + 62) 

One of the two integrations can be carried out, and leads to the following result [41] 



3 W 2 e-l™- m l«cos(n + 



mix 



G A A(n,m) = — dx : ' , (A.2) 



2tt Jo sin x 



COS 2 X 



where the function s(x) is defined through 



sinh s = SmX — cos 2 x. (A. 3) 

cosx 

The previous integral is still divergent, but provides a convergent integral representation for the 
subtracted Green function $AA{n,m) = GAA(n,m) — Gaa(0,0)- It yields the following values for 
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small n, m [42], 

QaaO-,0) = = ~, (A.4) 

^ AA (1, 2) = $aa(-1, 1) = 1 - — , (A.5) 

7T 

<S> AA (2, 0) = $^(2, 2) = -4 + (A.6) 

7T 

d>^(2, 3) = 2) = - - — , (A.7) 

2 7T 

^a(3,0) = $aa(3,3) = + — . (A.8) 

Let us now evaluate the asymptotic behaviour of <& AA (n, m) for large distances. We will do this 
calculation by using ideas from [41] and [13]; the analogous calculation for the square lattice has 
been done in [43]. 

The basic idea underlying these computations is that for large \n — m\, the exponential factor 
in (A. 2) contributes significantly only in the region where s is small, which is also where x is small. 
In this region, we may expand s(x) in powers of x, 



s(x) = VS[x+—x b + — x 9 + ...). (A.9) 

Therefore the main contribution of the integral comes from the part close to the origin and is given 
by the way the rest of the integrand behaves for small x. 

We start by splitting the integral giving & AA (n,m) into three pieces, 



3 rl 2 

<& AA (n,m) = — dx 

27T Jo 



vr/2 ^ e -| n-m|s cos ( n + m ) g _ 1 

sin x a/4 — cos 2 x 



3 r/ 2 c e -\ n - m \ s cos {n + m)x e -^\n-m\x cos ( n + m ) x 
2tt JO k SIM \/4 — cos 2 x V3x > 



\/4 — cos 2 x \/3a 
3 W 2 r e^^l™-" 1 ^ cos (n + m)x 1 
27ri ^ >/3x >/& 

3 W 2 r 1 1 



3 /-7T/2 I I 

+ — / dx{^= [. (A.10) 

27r Jo ^\/3z sin x v4 — cos 2 x > 

The second integral can be evaluated exactly, up to exponentially small terms, and turns out to 
give the dominant, logarithmic term, equal to — ^ (log r + 7 + log n) , where r 2 = n 2 + m 2 — nm and 
7 = 0.577216 is the Euler constant. The third integral is a constant which can also be computed 
exactly, and is equal to ^ log ^ . We obtain at this stage 



* 1 , ^ 
9 AA (n,m) = - — 

Z7T 



logr + 7 + - log 12 



3 d | e ~' n ~ m ' Scos ( n + m ) x e -v / 3|—k cos(n + m)x , (All) 

2vr JO SIM \/4 — cos 2 x V3x > 
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To evaluate the remaining integral, we use the expansion of s as a power series in x, and write 

e -\n- m \s = e -V3\n-m\xQ( x ) where Q ig expan ded as 



Q(x) = exp [— |n — m\(s — V3x)] = exp j — \/3| 
The subtracted Green function then becomes, with p 



n — m 



n 



z 5 z 9 

m| and q = n + m, 



(A.12) 
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(A.13) 



• smx — cos 2 x \/3x - 

By construction, the function in brackets is regular at x = 0, and may be expanded in powers of 
x. It is not difficult to see from (the polar coordinates have been introduced in Section 2) 



rn/2 

Jo 



dx e 



cos qx 



f 

Jo 



dx e 



cos qx 



v^p v^lcosc^- -jjsin^l 



io Jo 3p 2 + q 2 4r 

where we have neglected exponentially small terms, that the following estimate holds 

rvr/2 



I dx x k ~ 1 e~ y/3px cos qx ~ 0(r~ k ). 
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(A.14) 



(A.15) 



As a consequence, the integral in (A.13) has an expansion in inverse powers of r, for which the 
calculation of the r~ k terms requires the expansion of the function in brackets to order k — 1. 
Because Q(x) has coefficients which depend on p = 0(r), the order k — 1 means that we keep those 
terms p a x b such that b — a = k — 1. And since Q(x) is divided by sin x ~ x, Q(x) is to be expanded 
to order k. One easily checks that for fixed k, there is only a finite number of terms to consider. 
The rest of the calculations is straightforward. 

To order r~ 8 , the relevant expansion of Q(x) reads 



Q{X) = 1 " IT* "loT* + 675* 
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(A.16) 



from which we obtain the asymptotic expansion of the Green function 
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(A.17) 



We note that it is invariant under the symmetries of the lattice, generated by ip — > p + and 
ip — > I — if. When p = 0, corresponding to the line n = m or 99 = |, the above calculation breaks 
down. However this line is related by a symmetry of the lattice to <j> = it for which p = \n\ is not 
zero. The previous result is therefore valid for all p. 

As particular cases, we find the asymptotic behaviour on the line n = m (ip = |) 



$AA(m,m) = - 
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(A.18) 
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and on the line n = 2m (¥> = §), 

i 1 K 7 i 

(A.19) 
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For the intended calculations in the sandpile model, we also need to know the Green function for 
sites in the close neighborhood of a reference site. For this it is sufficient to compute <&AA(n + £, m + 
k) with m, n » k,£ as the other entries <&ab and <&ba may be obtained from them. To compute 
&AA(n + £,m + k), one may simply follow the above calculations in which one appropriately shifts 
n and m by £ and k respectively, and then expand the result in inverse powers of r. At order 4, we 
obtain, where r and 4> are the polar coordinates of the site (n, m) as before, 



/3 

$A4(™ + ^,m + fc) = - — 

Z7T 



log r + 7 + - log 12 



(21 - k) cos ip + v^/c 
47r r 



V3 (2£ 2 - 2k£ - k 2 ) cos 2ip + VSk(2£ - k) sin 2<p 
+ 87 7 2 

V3 (2£ 3 - 3k£ 2 - 3k 2 £ + 2k 3 ) cos 3ip + 3\/3A;£(^ - k) sin 3<p 
~ 12tt ^ 

y/3 4 cos 6^ + 15(2^ 4 - 4k£ 3 - 6k 2 £ 2 + 8k 3 £ - k A ) cos 4ip + 15\/3/c(4^ 3 - 6k£ 2 + A; 3 ) sin 4<^ 
+ 240vr ^ 

+ . . . (A.20) 
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